Mobility Restriction and the Spread of COVID-19 in the United States

Code

library(tidyverse)
library(lubridate)
library(urbnmapr)
library(zoo)
library(psych)
library(meta)
library(DT)
library(viridis)
library(gganimate)
library(gifski)

Overview

Mobility restriction is one of the primary policies implemented to mitigate the spread of the COVID-19 pandemic. However, its effectiveness remains a topic of debate. In this project, we aim to investigate the impact of mobility restrictions on the spread of COVID-19 in the U.S. The analysis considers various factors, such as government responses to COVID-19 and vaccination.

This assignment was conducted under the guidance of Dr. John H. Holmes, PhD, FACE, FACMI, Professor of Medical Informatics in Epidemiology.

My GitHub repo for this project can be found here

Introduction

Throughout the COVID-19 pandemic, human mobility restrictions were one of the main policies implemented in many countries with the aim of reducing the transmission of SARS-CoV-2.^1–3 Such restrictions include physical distancing and community containment measures to reduce public transport use, public gatherings, school closures, and encouraging working from home where possible. Prior experiences with the 2009 H1N1⁴ and Ebola⁵ provided evidence for the effectiveness of these interventions in reducing disease transmission. However, the effectiveness of imposing mobility restrictions as a policy for controlling COVID-19 outbreaks has been controversial. Recent studies found that travel restrictions were effective in the early stages of the outbreak but may be less useful when the disease is widespread.^6,7 Furthermore, these restrictions have also led to enormous economic losses.⁸ Estimates suggest that global GDP growth has fallen by as much as 10%, at least part of which can be attributed to these mobility restrictions.⁹ Hence, it is critical to quantify the effectiveness of decisions to apply large-scale mobility restrictions in limiting the spread of the pandemic.

This study was conducted to investigate the impact of mobility restrictions on the spread of COVID-19 in the U.S. In fact, the U.S. experienced significant challenges in managing the pandemic in the early stages due to the large population, diverse demographics, and different healthcare settings across states. The findings of this project might not only support the current response but also contribute valuable evidence to the scientific community, informing future disease control strategies. This evidence guides policymakers in future outbreaks, enabling timely interventions, especially when vaccination is not available.

Methods

Data sources

Table 1: Datasets Overview
Dataset	Measurement	Temporal coverage	Source
COVID cases	Daily	Jan 2020-present	usafacts.org
Google Mobility	Daily	Feb 2020-present	google.com
Vaccination and Government response	Daily	Jan 2020-present	ourworldindata.org

Table 1 summrizes the data sources used in this study

COVID-19 cases and deaths

We obtained the data for COVID-19 cases and deaths from the USAFacts official website (https://usafacts.org/visualizations/coronavirus-covid-19-spread-map). USAFacts is a non-profit civic initiative that aims to provide a data-driven portrait of the American population, U.S. government finances, and the government’s impact on society. For COVID-19, USAFacts offers real-time pandemic data from all 50 states and the capital city.¹⁰

Google Mobility Dataset

For mobility measures, we used the Google mobility dataset. The data contain information about the daily amount of visitors to a specific place, including (1) groceries and pharmacies, (2) transit stations, (3) retail and recreation venues, (4) workplaces, (5) parks, and (6) residential areas. The measurements were based on mobile device-based global positioning system (GPS).¹¹ The data were measured from February 15, 2020 to present.

Vaccination and Government responses

We collected data on vaccination and Government responses from Our World in Data (OWID). The OWID is managed by a non-profit organization, providing rich datasets for the COVID-19 pandemic over the globe, including cases, deaths, vaccination, policies, population characteristics, etc.¹²

Code

#---------- COVID cases and deaths
df_case <- read_csv("Data/covid_confirmed_usafacts.csv")
df_death <- read_csv("Data/covid_deaths_usafacts.csv")

#----- Case: aggregate at state level
df_case <- df_case |>
  rename(county_name = `County Name`) |>
  gather(-c(countyFIPS, county_name, State, StateFIPS), 
            key = "date", value = "new_cases") |>
  mutate(date = ymd(date)) |> 
  group_by(State, date) |>
  summarise(new_cases_cum = sum(new_cases, na.rm = T)) |>
  ungroup() |>
  group_by(State) |>
  mutate(new_cases_cum_lag1 = lag(new_cases_cum, 1),
         new_cases = new_cases_cum - new_cases_cum_lag1,
         new_cases = ifelse(new_cases < 0, 0, new_cases)) |>
  ungroup()|>
  select(-new_cases_cum_lag1)

#----- Deaths: aggregate at state level
df_death <- df_death |>
  rename(county_name = `County Name`) |>
  gather(-c(countyFIPS, county_name, State, StateFIPS), 
         key = "date", value = "new_deaths") |>
  mutate(date = ymd(date))|> 
  group_by(State, date) |>
  summarise(new_deaths_cum = sum(new_deaths, na.rm = T)) |>
  ungroup()|>
  group_by(State) |>
  mutate(new_deaths_cum_lag1 = lag(new_deaths_cum, 1),
         new_deaths = new_deaths_cum - new_deaths_cum_lag1,
         new_deaths = ifelse(new_deaths < 0, 0, new_deaths)) |>
  ungroup() |>
  select(-new_deaths_cum_lag1)

#----- Merge new case, death, and population datasets
df_oc <- df_case |> left_join(df_death, by = c("State", "date")) 

# Get counties shapefile: to get full state name from abbriviation
st_name_full <- urbnmapr::states |>
  group_by(state_abbv, state_name) |>
  slice(1) |>
  select(state_abbv, state_name)

# Create full name for state
df_oc <- df_oc |> rename(state_abbv = State) |>
  left_join(st_name_full, by = "state_abbv")

Code

#---------- Google Mobility data
df_mob <- read_csv("Data/Global_Mobility_Report.csv")

df_mob <- df_mob |> filter(country_region == "United States") |> 
  filter(sub_region_1 != "") |>
  rename(
    state_name = sub_region_1,
    grocery_pharm = grocery_and_pharmacy_percent_change_from_baseline,
    retail_recreation = retail_and_recreation_percent_change_from_baseline,
    park = parks_percent_change_from_baseline,
    transit = transit_stations_percent_change_from_baseline,
    workplace = workplaces_percent_change_from_baseline,
    residential = residential_percent_change_from_baseline,
  ) |>
  select(state_name, date, grocery_pharm, retail_recreation, park,
         transit, workplace, residential)

# Aggregated at state level
df_mob <- df_mob |> group_by(state_name, date) |>
  summarise(grocery_pharm = mean(grocery_pharm, na.rm = T)/100,
            retail_recreation = mean(retail_recreation, na.rm = T)/100,
            retail_recreation = mean(retail_recreation, na.rm = T)/100,
            park = mean(park, na.rm = T)/100,
            transit = mean(transit, na.rm = T)/100,
            workplace = mean(workplace, na.rm = T)/100,
            residential = mean(residential, na.rm = T)/100)

Code

#---------- Vaccination
# Data from Our World in Data
df_owid <- read_csv("Data/owid-covid-data.csv") |>
  filter(location == "United States")

df_owid <- df_owid |>
  select(date, reproduction_rate, new_tests_smoothed_per_thousand, tests_per_case,
         total_vaccinations_per_hundred, people_fully_vaccinated_per_hundred,
         stringency_index) |>
  rename(
    test_thousand = new_tests_smoothed_per_thousand,
    vacc_any = total_vaccinations_per_hundred,
    vacc_fully = people_fully_vaccinated_per_hundred,
  )

Code

#---------- Merge all data
# setdiff(df_oc$state_name, df_mob$state_name)

df <- df_oc |> left_join(df_mob, by = c("state_name", "date")) |> 
  left_join(df_owid, by = "date")

# Exclude Omicron variant
df <- df |>
  filter(date >= ymd("2020-02-15") & date < ymd("2021-12-01"))

# Replace missing values (testing, vaccine were not available)
df <- df |>
  mutate(
    test_thousand = ifelse(is.na(test_thousand), 0, test_thousand),
    tests_per_case = ifelse(is.na(tests_per_case), 0, tests_per_case),
    vacc_any = ifelse(is.na(vacc_any), 0, vacc_any),
    vacc_fully = ifelse(is.na(vacc_fully), 0, vacc_fully)
  )

dim(df)

[1] 33405    19

Code

df[1:200, ] |> datatable()

Variables

Outcomes

The main outcome of this study was the growth rate (GR) of cases. The GR indicates how fast the spread of COVID-19 is. In this study, we defined GR of cases for specific state ith and date tth as the logarithmic rate of change for the new cases (C) in the preceding three days relative to the logarithmic rate of change for the new cases in the preceding seven days.

\[GR_i^t = \frac{log(\sum_{t-2}^{t}\frac{C_i^t}{3})}{log(\sum_{t-6}^{t}\frac{C_i^t}{7})}\]

Independent variables

The main exposure in this study was mobility, derived from the original values of the Google mobility dataset. These mobility values reflect the median change (in percentage) in the number of visitors to specific categories of locations compared to the reference period (between January 3 and February 6, 2020, before the declaration of COVID-19 as a global pandemic).¹¹ The data were aggregated at the state level.

We defined the main mobility variable by average mobility of six venues (groceries and pharmacies, transit stations, retail and recreation venues, workplaces, parks, and residential areas)

Covariates

The covariates included

Fully vaccination, which was defined as the percentage of the population having fully two required doses of COVID-19 vaccination.
Stringency index, which was a composite measure of nine of the response metrics: school closures; workplace closures; cancellation of public events; restrictions on public gatherings; closures of public transport; stay-at-home requirements; public information campaigns; restrictions on internal movements; and international travel controls.¹² The values of the stringency index were calculated as the mean scores of the nine metrics, with the range of 0 to 100. A higher value of the stringency index indicates a stricter response.¹²

Code

#---------- Outcome
df <- df |> 
  group_by(state_name) |>
  mutate(
    GR_case_num = lag(rollapply(new_cases, 3, mean, fill = NA), 1),
    GR_case_den = lag(rollapply(new_cases, 7, mean, fill = NA), 3),
    GR_case = log(GR_case_num)/log(GR_case_den)
  ) |> ungroup() 


df <- df |>
  mutate(GR_case = ifelse(is.nan(GR_case), 0, GR_case)) |>
  mutate(GR_case = ifelse(is.infinite(GR_case), NA, GR_case))

Code

#---------- Composite mobility
df <- df |> mutate(
  mobility = (grocery_pharm + transit + workplace + 
                     retail_recreation + park + residential)/6
)

df[1:200, ] |> datatable()

Statistical analysis

Descriptive statistics

For descriptive statistics, we created the animation choropleth maps to show the GR and six categories of mobility for all 50 states and the capital city over time. The scatter plots were used to visualize the initial sign of the correlation between the overall growth rate and outcome across all states oever time.

Effects of mobility restrictions on spread of COVID-19

To evaluate the effects of mobility restrictions on the spread of COVID-19, we first restricted the data on the date before omicron variant was detected (i.e., November 30, 2021).

Multiple linear regressions were performed to estimate the relationship of interest. The models were fitted separately for each state. The random effect meta-analysis was used to pool the estimates across 50 states and the capital city. The \(I^2\) statistic was used to evaluate the heterogeneity in the estimates across states.

As mobility required delayed time to affect the GR of cases of COVID-19, we applied the lag effect of mobility. Based on the previous study, the lag of 14 days was considered as the optimal lag for mobility on GR of cases.⁷ Therefore, we selected a lag of 14 days in our main analysis, which implied that mobility policy today would affect GR of cases in the following 14 days.

In the main analysis, the composite mobility was obtained by taking the average of six mobility indicators. In the subsequent analyses, each of the six mobility indicators was analyzed.

All the models were adjusted for full vaccination and stringency index. All these confounders were applied with the lag effects of 14 days.

Stratified analysis

We assume that the effects of mobility restrictions on the spread of COVID-19 are different across the periods of COVID-19 pandemic. Therefore, we conducted the analysis stratified by two periods: (1) from February 15, 2020 (when the first mobility measure was available) to the date of 10% of the population received fully vaccination (before vaccine period), (2) after the date of 10% of the population received fully vaccination (after vaccine period).

Sensitivity analysis

We did sensitivity analysis to evaluate the relationship of interest under different lag effects and different types of mobility. Therefore, we repeated the analyses with different lag, from 1 day to 21 days prior to the GR, for each type mobility.

Results

Descriptive statistics

Animation maps for outcome and mobility over time

Code

# Summarized by month for animation maps
df2 <- df |>
  mutate(month_year = ym(paste0(year(date), "-", month(date))))


df2 <- df2 |> group_by(state_abbv, state_name, month_year) |>
  summarise(
    GR_case = mean(GR_case, na.rm = T),
    mobility = mean(mobility, na.rm = T),
    grocery_pharm = mean(grocery_pharm, na.rm = T),
    transit = mean(transit, na.rm = T),
    workplace = mean(workplace, na.rm = T),
    park = mean(park, na.rm = T),
    retail_recreation = mean(retail_recreation, na.rm = T),
    residential = mean(residential, na.rm = T),
  ) |>
  ungroup()

# Merge data with shape file for US
df2_sh <- df2 |>
  left_join(urbnmapr::states, by = c("state_abbv", "state_name"))

# Set up theme
my_theme <- function() {
  theme_minimal() +                                  
  theme(axis.line = element_blank(),                 
        axis.text = element_blank(),                 
        axis.title = element_blank(),
        panel.grid = element_line(color = "white"),  
        legend.key.size = unit(0.4, "cm"),          
        legend.text = element_text(size = 16),       
        legend.title = element_text(size = 16),
        plot.title = element_text(size = 20),
        legend.position = "bottom",
        strip.text = element_text(face = "bold", size = 16),
        strip.background = element_rect(fill = "white", color = NA))      
}

# Growth rate
#--------------------------------------------------------------
GR_trans_fig <- df2_sh |> filter(GR_case > 0) |>
  ggplot() + 
  geom_polygon(aes(long, lat, group = group, fill = GR_case), 
               color = "white", linewidth = 0.02) +
  coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
  scale_fill_viridis(option = "D",
                     direction = -1, 
                     name = "Growth rate",
                     guide = guide_colorbar(
                     direction = "horizontal",
                     barheight = unit(2, units = "mm"),
                     barwidth = unit(100, units = "mm"),
                     draw.ulim = FALSE,
                     title.position = "top",
                     title.hjust = 0.5,
                     title.vjust = 0.5)) +
  my_theme() +
  labs(
    title = "Growth rate of cases on {frame_time}"
  ) + 
  transition_time(month_year)

# Mobility
#--------------------------------------------------------------
mob_trans_fig <- df2_sh |> filter(GR_case > 0) |>
  ggplot() + 
  geom_polygon(aes(long, lat, group = group, fill = mobility*100), 
               color = "white", linewidth = 0.02) +
  coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
  scale_fill_viridis(option = "C",
                     direction = -1, 
                     name = "Mobility (%)",
                     guide = guide_colorbar(
                     direction = "horizontal",
                     barheight = unit(2, units = "mm"),
                     barwidth = unit(100, units = "mm"),
                     draw.ulim = FALSE,
                     title.position = "top",
                     title.hjust = 0.5,
                     title.vjust = 0.5)) +
  my_theme() +
  labs(
    title = "Mobility on {frame_time}"
  ) + 
  transition_time(month_year)

Code

animate(GR_trans_fig)

Figure 1: Growth rate of cases across states over time

Code

animate(mob_trans_fig)

Figure 2: Overall mobility across states over time

Figure 1 and Figure 2 describe the time series changes of GR of cases and mobility across the 50 states and the capital city. Overall, there is a noticeable increase in the GR over time, although a slight drop was observed in a short period at the beginning of 2021.

The overall mobility varied over time, with the most restricted period (i.e., the period showing the smallest values of mobility) occurring between September 2021 and January 2022.

Animation maps for six types of mobility

Code

df2_sh_long <- df2_sh |> select(grocery_pharm, transit, workplace, park, 
                                retail_recreation, residential, long, lat, 
                                group, month_year) |>
  gather(-c(long, lat, group, month_year), key = "mob_type", value = mobility) |>
  mutate(mob_type2 = case_when(mob_type == "grocery_pharm" ~ "Groceries & pharmacies",
                               mob_type == "transit" ~ "Transit stations",
                               mob_type == "workplace" ~ "Workplaces",
                               mob_type == "retail_recreation" ~ "Retail and recreation",
                               mob_type == "park" ~ "Parks",
                               mob_type == "residential" ~ "Residential areas")) 

mob_all_trans_fig <- df2_sh_long |>
  ggplot() + 
  geom_polygon(aes(long, lat, group = group, fill = mobility*100), 
               color = "white", linewidth = 0.02) +
  coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
  scale_fill_viridis(option = "C",
                     direction = -1, 
                     name = "Mobility (%)",
                     guide = guide_colorbar(
                     direction = "horizontal",
                     barheight = unit(2, units = "mm"),
                     barwidth = unit(100, units = "mm"),
                     draw.ulim = FALSE,
                     title.position = "top",
                     title.hjust = 0.5,
                     title.vjust = 0.5)) +
  facet_wrap(~mob_type2, ncol = 3) +
  my_theme() +
  labs(
    title = "Mobility on {frame_time}"
  ) +
  transition_time(month_year)

animate(mob_all_trans_fig, width = 1200, height = 700)

Figure 3: Six types of mobility across states over time

Among the six types of mobility, movement to parks exhibited the most variation, while mobility to retail and recreation venues showed the least variation over time (see Figure 3).

For a more detailed depiction of the variations in average values across states concerning GR, overall mobility, and the six types of mobility over time, please see Figure 4 and Figure 5.

Code

df <- df |> 
  mutate(period = ifelse(vacc_fully < 10, 0, 1),
         period = factor(period, labels = c("Before vaccine", "After vaccine")))

df_long <- df |> select(grocery_pharm, transit, workplace, park,
                            retail_recreation, residential, date) |>
  gather(-date, key = "mob_type", value = mobility) |>
  mutate(
    mob_type2 = case_when(mob_type == "grocery_pharm" ~ "Groceries & pharmacies",
                          mob_type == "transit" ~ "Transit stations",
                          mob_type == "workplace" ~ "Workplaces",
                          mob_type == "retail_recreation" ~ "Retail and recreation",
                          mob_type == "park" ~ "Parks",
                          mob_type == "residential" ~ "Residential areas"),
    mob_type2 = factor(mob_type2, 
                      levels = c("Groceries & pharmacies", 
                                 "Transit stations", "Workplaces", 
                                 "Retail and recreation", 
                                 "Parks", "Residential areas"))) 


fig_des_GR <- df |>
  ggplot(aes(x = date, y = GR_case)) +
  geom_jitter(alpha = 0.2, color = "#016c59", shape = 21, size = 0.5) +
  geom_smooth(se = F, color = "#feb24c") +
  labs(x = "Date", y = "Growth rate of cases", 
       title = "Growth rate of cases") +
  theme_minimal() + 
  theme(
    legend.text = element_text(size = 16),       
    legend.title = element_text(size = 16),
    plot.title = element_text(size = 20)
  ) +
  coord_cartesian(ylim = c(-2, 2)) 

fig_des_mob_overall <- df |>
  ggplot(aes(x = date, y = mobility*100)) +
  geom_jitter(alpha = 0.2, color = "#e7298a", shape = 21, size = 0.5) +
  geom_smooth(se = F, color = "#feb24c") +
  labs(x = "Date", y = "Mobility (%)", title = "Overall mobility") +
  theme_minimal() + 
  theme(
    legend.text = element_text(size = 16),       
    legend.title = element_text(size = 16),
    plot.title = element_text(size = 20)
  ) +
  coord_cartesian(ylim = c(-50, 150)) 

cowplot::plot_grid(fig_des_GR, fig_des_mob_overall, ncol = 1, labels = "AUTO")

Figure 4: Average mobility and growth rate of cases over time

Code

fig_des_mob <- df_long |>
  ggplot(aes(x = date, y = mobility*100)) +
  geom_smooth(aes(color = mob_type2, fill = mob_type2)) +
  scale_color_brewer(palette = "Set1") +
  scale_fill_brewer(palette = "Set1") +
  labs(x = "Date", y = "Mobility (%)",
       title = "Six types of mobility",
       color = NULL,
       fill = NULL) +
  theme_minimal() +
  theme(legend.text = element_text(size = 14), 
        plot.title = element_text(size = 20),
        panel.grid.minor = element_blank(),
        legend.position = "bottom") +
  coord_cartesian(ylim = c(-50, 150)) 

fig_des_mob

Figure 5: Average (six types) mobility over time

Effects of mobility restrictions on spread of COVID-19

Main analysis

To reduce the repeated work across states (50 + 1), different lag effects of mobility (ranging from 7 to 21 days), multiple types of mobility (overall and six specific types), and stratification into two periods (pre and post-vaccine), I developed a function to estimate the effects of interest. This function was used for subsequent analyses.

Code

# Obtain state name for loop
st_name <- df$state_name |> unique()
lag_name <- paste0(rep(c("beta", "se"), 15), rep(c(7:21), each = 2))

# Define the function
make_state_est <- function(var, ...) {
  
  # Define matrices to store the results
  mat_earl <- matrix(ncol = 30, nrow = length(st_name))
  mat_late <- matrix(ncol = 30, nrow = length(st_name))
  
  # 2 nested loop for `states` and different `lags`
  for (i in seq(length(st_name))) { 
    # Data for each state, here I convert `tible` to normal `data.frame` 
    # to run easier inside a loop
    df_sub <- as.data.frame(subset(df, state_name == st_name[i]))
    df_sub$vacc_fully_lag <- lag(df_sub$vacc_fully, 14)
    df_sub$stringency_index_lag <- lag(df_sub$stringency_index, 14)
    # Apply different lags 7-21 days
    for (j in 1:15) {
      df_sub$mobility_lag <- lag(df_sub[, var], j+6)
      
      # Fit model 
      # Early period
      m1 <- lm(GR_case ~ mobility_lag + vacc_fully_lag + stringency_index_lag, 
               data = subset(df_sub, vacc_fully_lag < 10))
      
      # Late period
      m2 <- lm(GR_case ~ mobility_lag + vacc_fully_lag + stringency_index_lag, 
               data = subset(df_sub, vacc_fully_lag >= 10))
      
      m1_summ <- summary(m1)
      m2_summ <- summary(m2)
      
      # Store results to 2 matrices
      mat_earl[i, j*2-1] <- m1_summ$coefficients[2, 1]
      mat_earl[i, j*2] <- m1_summ$coefficients[2, 2]
      
      mat_late[i, j*2-1] <- m2_summ$coefficients[2, 1]
      mat_late[i, j*2] <- m2_summ$coefficients[2, 2]
      
    }
  }
  # Convert to dataframes
  mat_earl <- as.data.frame(mat_earl)
  mat_late <- as.data.frame(mat_late)
  
  names(mat_earl) <- lag_name
  names(mat_late) <- lag_name
  
  mat_earl$state_name <- st_name
  mat_late$state_name <- st_name
  
  return(list(early = mat_earl, late = mat_late))
}

Here I apply the function to obtain two estimates across all states and all 15 lag effects for early and late periods

Code

# Apply the function
mobi_df_list <- make_state_est(var = "mobility")

# Store two datasets for later analysis
ovarall_earl <- mobi_df_list$early
ovarall_late <- mobi_df_list$late

The meta-analysis showing the effects of overall mobility on the GR of COVID-19 cases before and after the vaccination periods is presented in Figure 6 and Figure 7.

During the pre-vaccination period, mobility demonstrated a positive association with the GR of COVID-19 (\(\beta\) = 0.11, 95% CI: 0.07 to 0.15). This suggests that restricting mobility could potentially reduce the spread of COVID-19 (see Figure 6).

Conversely, during the post-vaccination period, this relationship reversed, showing \(\beta\) = -0.07 (95% CI: -0.10 to -0.03). This indicates that mobility restrictions during this phase might not be as effective and could potentially increase the spread of COVID-19 (see Figure 7).

Note that, the results were obtained with 14-day lag effect of mobility.

Code

# Conduct meta-analysis to pool the effect of all states
mobility_earl_meta <- metagen(TE = beta14,
                          seTE = se14,
                          studlab = state_name,
                          method.tau = "DL",
                          common = FALSE,
                          data = ovarall_earl,
                          title = "Effect of mobility on growth rate, before vaccine period")
# Forest plot
forest(mobility_earl_meta)

Figure 6: Relationship between mobility and growth rate of cases, with 14-day lag effect and before vaccine period

Code

mobility_late_meta <- metagen(TE = beta14,
                          seTE = se14,
                          studlab = state_name,
                          method.tau = "DL",
                          common = FALSE,
                          data = ovarall_late,
                          title = "Effect of mobility on growth rate, after vaccine")
# Forest plot
forest(mobility_late_meta)

Figure 7: Relationship between overall mobility and growth rate of cases, with 14-day lag effect and after vaccine period

Sensitivity analysis

Sensitivity analsysis for different lag effects of mobility (7 to 21 days)

As I proposed sensitivity analyses for all types of mobility, the function was created for this purpose.

Code

sens_analysis <- function(data_early, data_late, title = NA, limit) {
  
  # Define matrices to store the results
  early_sen_mat <- matrix(ncol = 2, nrow = 15)
  late_sen_mat <- matrix(ncol = 2, nrow = 15)
  
  # loop for different lags
  for (i in 1:15) {
    meta_early_i <- metagen(TE = data_early[, i*2-1],
                            seTE = data_early[, i*2],
                            method.tau = "DL",
                            common = FALSE)
    meta_late_i <- metagen(TE = data_late[, i*2-1],
                           seTE = data_late[, i*2],
                           method.tau = "DL",
                           common = FALSE)
    
    # Store result to matrices
    early_sen_mat[i, 1] <- meta_early_i$TE.random
    early_sen_mat[i, 2] <- meta_early_i$seTE.random
    
    late_sen_mat[i, 1] <- meta_late_i$TE.random
    late_sen_mat[i, 2] <- meta_late_i$seTE.random
  }
  
  early_sen_mat <- as.data.frame(early_sen_mat)
  late_sen_mat <- as.data.frame(late_sen_mat)
  
  names(early_sen_mat) <- c("pooled_beta", "se")
  early_sen_mat$lag <- 7:21
  
  names(late_sen_mat) <- c("pooled_beta", "se")
  late_sen_mat$lag <- 7:21
  
  mobility_sen_mat <- rbind(early_sen_mat, late_sen_mat) |>
    mutate(low = pooled_beta - 1.96*se,
           high = pooled_beta + 1.96*se)
  
  mobility_sen_mat$period <- rep(c("Before vaccine", "After vaccine"), each = 15)
  
  
  # Sensitivity figure
  dodge <- position_dodge(width = 0.5)
  
  fig_sen <- ggplot(aes(as.factor(lag), y = pooled_beta, 
                        ymin = low, ymax = high,
                        color = period), data = mobility_sen_mat) +
    geom_errorbar(width = 0.4, linewidth = 1, position = dodge) +
    geom_point(size = 2, shape = 21, fill="white", position = dodge) +
    geom_hline(yintercept = 0, linetype = 2, col = "gray10", size = 1) +
    theme_minimal() +
    scale_color_brewer(palette = "Set1") +
    coord_cartesian(ylim = limit) +
    theme_bw() + 
    theme(legend.position = "top") +
    labs(
      x = "Lag of mobility (day)",
      y = "Pooled estimate of all states",
      title = title,
      color = NULL
    )
  
  return(fig_sen)
}

Apply the function to evaluate effects of different lags of mobility on growth rate of cases

Code

sens_analysis(ovarall_earl, ovarall_late, 
              title = "Effect of mobility (overall) on growth rate", 
              limit = c(-0.3, 0.3))

Figure 8: Sensitivity analysis for the relationship between overall mobility and growth rate of cases

The sensitivity analysis examining various lag effects of overall mobility on the GR of COVID-19 is depicted in Figure 8.

During the pre-vaccine period, consistent effects were observed across different lag periods, consistently indicating a positive impact of mobility restriction on the GR. Notably, the optimal lag effects were identified between 12 to 15 days of mobility. This suggests that mobility policy implementations could yield the most significant impact on GR in the subsequent 12 to 15 days.

In contrast, during the later period when vaccines were available, the effects showed inconsistency. This suggests that the mobility restriction policy might not be as effective during this phase (Figure 8).

Sensitivity analsysis for different types of mobility

In this section, I applied the two functions (i.e., make_state_est and sens_analysis) defined in the Section 5.2.1. Beside, I created a loop to conduct all sensitivity analyses for 6 types of mobility simultaneously.

Code

mob_label <- c("groceries & pharmacies", "transit stations", "workplaces",
               "retail and recreation", "parks", "residential areas")
mob_name <- c("grocery_pharm", "transit", "workplace", 
              "retail_recreation", "park", "residential")
limit_range <- c(rep(c(-0.3, 0.3), 3), c(-15, 15), c(-0.05, 0.05), c(-1, 1))

# For all 6 types of mobility
for (i in 1:6) {
  data_list <- make_state_est(var = mob_name[i])
  zzz <- paste0("fig_", mob_name[i])
  eval(call("<-", as.name(zzz), 
            sens_analysis(data_list$early, data_list$late, 
              title = paste0("Effect of mobility for ", mob_label[i], " on growth rate"), 
              limit = c(limit_range[i*2-1], limit_range[i*2]))
              ))
}

# Plot all figures
cowplot::plot_grid(fig_grocery_pharm, fig_transit, fig_workplace,
                   fig_retail_recreation, fig_park, fig_residential,
                   labels = "AUTO", ncol = 2)

Figure 9: Sensitivity analysis for the relationship between six types of mobility and growth rate of cases

The similar patterns were observed for groceries & pharmacies, transit stations, workplaces, retail and recreation, and parks (see Figure 9).

However, the effects of mobility in residential areas were contradictory. The estimates for the relationship between this mobility category and GR of COVID-19 were primarily negative during most lag effects in the pre-vaccine period.

This can be attributed to the nature of this mobility category. If mobility restrictions were applied, reducing movement to public spaces, transit, and workplaces, the population tended to stay within residential areas, resulting in effects opposite to those observed in the previous five categories.

Discussion and conclusion

In this project, we investigated the impact of mobility restrictions on the spread of COVID-19 in the U.S.

We found that, during the pre-vaccination period, mobility was positively associated with GR of COVID-19, suggesting that mobility restrictions could potentially reduce the spread of COVID-19. However, in the post-vaccination period, this relationship reversed.

This reversal might be attributed to the post-vaccination period leading individuals to be less diligent in practicing self-protection methods, such as physical distancing and wearing face masks. Consequently diminishing the impact of mobility restrictions.

The sensitivity analyses demonstrating the effects during the pre-vaccine period were consistent, but not for the post-vaccine period.

We found the optimal lag effects of mobility were between 12 to 15 days during pre-vaccine period.

Analyses for five types of mobility (groceries & pharmacies, transit stations, workplaces, retail and recreation, and parks) yielded similar findings to overall mobility. However, mobility in residential areas exhibited contrasting estimates, which might be attributed to its distinctive nature.

Limitations

The results of this study should be interpreted in the context of potential limitations. Firstly, the spread of COVID-19 is influenced by various self-protection behaviors (like physical distancing, mask-wearing, and hand hygiene), which we lacked data on, and thus we couldn’t control for those factors in our model. Secondly, our utilization of Google mobility data was based solely on individuals who shared their locations, potentially limiting the representation of the entire US population’s mobility patterns. Lastly, the seasonal variability of mobility, for example, the difference between summer and winter, was not fully captured as the mobility data were only available from early 2020. Ideally, using reference periods from the same weeks in previous years as baselines would be more precise than relying solely on the period between January 3 and February 6, 2020.

Recommendation

Based on our findings, we recommend implementing suitable mobility restriction policies during the early period when vaccines are unavailable. Further research is necessary to explore the impact of mobility restrictions during the later period, requiring more comprehensive individual-level data, particularly regarding self-protection behaviors.

References

Flaxman, S. et al. Estimating the effects of non-pharmaceutical interventions on COVID-19 in europe. Nature 584, 257–261 (2020).

The Lancet. India under COVID-19 lockdown. Lancet (London, England) 395, 1315 (2020).

Signorelli, C., Scognamiglio, T. & Odone, A. COVID-19 in italy: Impact of containment measures and prevalence estimates of infection in the general population. Acta Bio Medica: Atenei Parmensis 91, 175 (2020).

Chowell, G. et al. Characterizing the epidemiology of the 2009 influenza a/H1N1 pandemic in mexico. PLoS medicine 8, e1000436 (2011).

Peak, C. M. et al. Population mobility reductions associated with travel restrictions during the ebola epidemic in sierra leone: Use of mobile phone data. International journal of epidemiology 47, 1562–1570 (2018).

Kraemer, M. U. et al. The effect of human mobility and control measures on the COVID-19 epidemic in china. Science 368, 493–497 (2020).

Oh, J. et al. Mobility restrictions were associated with reductions in COVID-19 incidence early in the pandemic: Evidence from a real-time evaluation in 34 countries. Scientific reports 11, 13717 (2021).

Han, E. et al. Lessons learnt from easing COVID-19 restrictions: An analysis of countries and regions in asia pacific and europe. The Lancet 396, 1525–1534 (2020).

Hatzius, J. et al. The sudden stop: A deeper trough, a bigger rebound. Goldman Sachs Economics Research 31, (2020).

10.

USAFacts. DUS COVID-19 cases and deaths by state, https://usafacts.org/visualizations/coronavirus-covid-19-spread-map, Accessed: 11-10-2023.

11.

Google. COVID-19 Community Mobility Reports, https://www.google.com/covid19/mobility/, Accessed: 11-10-2023.

12.

Our World in Data. Coronavirus Pandemic (COVID-19), https://ourworldindata.org/coronavirus, Accessed: 11-10-2023.